/*==============================================================================
FILE NAME: Figure_3.do
CREATED: 12 June 2025
==============================================================================*/

**Figure 3

/* Set directory if working independently through code
if c(username)=="" { //insert username
	global rootdir "" // insert root path
	global processed_data "$rootdir/processed_data" 
	global figures "$rootdir/output/figures"  // Define global paths for replication package
} 
*/

//Panel A
//any air inv
set scheme modern
use "$processed_data/Air_Panel.dta", clear

* Drop RNs that never had air
drop if never_air_inv == 1

* Create time variables
egen t = group(year month)
xtset RN_id t

* Generate forward difference variables (0-12 months after incident)

forv h = 0/12 {
	gen p_air_inv_`h' = f`h'.p_air_inv - l1.p_air_inv
}

* Generate backward differences (2-12 months before incident)
forv h = 2/12 {
	gen p_air_inv_neg`h' = l`h'.p_air_inv - l1.p_air_inv
}

* RN-year fixed effects
egen RN_year = group(RN year)

* Keep relevant variables
keep RN_id t RN_year p_air_inv* p_air_incident 

* Initialize storage variables for regression results
cap drop b u d se Years Zero
gen Years = _n-13 if _n<=25
gen Zero = 0 if _n <=25 
gen b = 0
gen se = 0
gen u = 0 
gen d = 0

* Run regressions for likelihood of an investigation before and after citizen complaints
foreach h in 0 1 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_inv_`h'  p_air_incident, absorb(RN_year t) cluster(RN_id )
unique t if e(sample)
replace b = _b[p_air_incident] if Years == `h'
replace se = _se[p_air_incident] if Years == `h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == `h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == `h'
}
foreach h in 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_inv_neg`h'  p_air_incident, absorb(RN_year t) cluster(RN_id)
unique t if e(sample)
replace b = _b[p_air_incident] if Years == -`h'
replace se = _se[p_air_incident] if Years == -`h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == -`h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == -`h'
}

keep if Years != .
keep b u d se Years Zero

* Create point estimates output
export delimited "$point_estimates/Point_Estimates_Figure_3_Panel_A.csv", replace

* Plot event study graph (Panel A)
graph set window fontface "Times New Roman"
twoway(rarea u d Years, col(gs10) fint(inten30) lwidth(0) lpattern(solid)) ///
    (line b Years, lcolor(gs3) lpattern(solid) lwidth(medium)) ///
    (line Zero Years, lcolor(gs8)), ///
    xlabel(-12(1)12, nogrid labsize(vlarge)) ///
    ylabel(0(0.1)0.6, labsize(vlarge)) ///
    legend(off) ///
    ytitle("{&Delta} P(Investigation)", size(vlarge)) ///
    xtitle("Month", size(vlarge)) ///
    graphregion(color(white)) ///
    plotregion(color(white)) ///
    xsize(8.6)
graph export "$figures/Figure_3_Panel_A.pdf", replace

* Panel B
use "$processed_data/Air_Panel_with_referred.dta", clear
keep RN_id t RN_year p_air_inv* p_referred p_air_incident

* Initialize variables
cap drop b u d se b_ref se_ref u_ref d_ref Years Zero
gen Years = _n-13 if _n<=25
gen Zero = 0 if _n <=25
gen b = 0
gen se = 0
gen u = 0 
gen d = 0
gen b_ref = 0
gen se_ref = 0
gen u_ref = 0 
gen d_ref = 0


* Run regressions for likeihood of investigation before and afer complaints recieved that are outside of TCEQs jurisdiction
foreach h in 0 1 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_inv_`h' p_referred p_air_incident, absorb(RN_year t) cluster(RN_id )
unique t if e(sample)
replace b = _b[p_air_incident] if Years == `h'
replace se = _se[p_air_incident] if Years == `h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == `h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == `h'
replace b_ref = _b[p_referred] if Years == `h'
replace se_ref = _se[p_referred] if Years == `h'
replace u_ref = (_b[p_referred] + 1.96*_se[p_referred]) if Years == `h'
replace d_ref = (_b[p_referred] - 1.96*_se[p_referred]) if Years == `h'
}
foreach h in 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_inv_neg`h' p_referred p_air_incident, absorb(RN_year t) cluster(RN_id)
unique t if e(sample)
replace b = _b[p_air_incident] if Years == -`h'
replace se = _se[p_air_incident] if Years == -`h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == -`h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == -`h'
replace b_ref = _b[p_referred] if Years == -`h'
replace se_ref = _se[p_referred] if Years == -`h'
replace u_ref = (_b[p_referred] + 1.96*_se[p_referred]) if Years == -`h'
replace d_ref = (_b[p_referred] - 1.96*_se[p_referred]) if Years == -`h'
}

* Save point estimates output
keep if Years != .
keep b u d se b_ref se_ref u_ref d_ref Years Zero
export delimited "$point_estimates/Point_Estimates_Figure_3_Panel_B.csv", replace

* Plot event study graph (Panel B)
graph set window fontface "Times New Roman"
twoway(rarea u_ref d_ref Years, col(gs10) fint(inten30) lwidth(0) lpattern(solid)) ///
    (line b_ref Years, lcolor(gs3) lpattern(solid) lwidth(medium)) ///
    (line Zero Years, lcolor(gs8)), ///
    xlabel(-12(1)12, nogrid labsize(vlarge)) ///
    ylabel(0(0.1)0.6, labsize(vlarge)) ///
    legend(off) ///
    ytitle("{&Delta} P(Investigation)", size(vlarge)) ///
    xtitle("Month", size(vlarge)) ///
    graphregion(color(white)) ///
    plotregion(color(white)) ///
    xsize(8.6)
graph export "$figures/Figure_3_Panel_B.pdf", replace


// Panel C
use "$processed_data/Air_Panel.dta", clear

drop if never_air_inv == 1

egen t = group(year month)
xtset RN_id t

* Generate difference variables
forv h = 0/12 {
	gen p_air_complaint_inv_`h' = f`h'.p_air_complaint_inv - l1.p_air_complaint_inv
}

forv h = 2/12 {
	gen p_air_complaint_inv_neg`h' = l`h'.p_air_complaint_inv - l1.p_air_complaint_inv
}

egen RN_year = group(RN year)

* Keep particular variables
keep RN_id t RN_year p_air_complaint_inv* p_air_incident

* Initialize storage
cap drop b u d se Years Zero
gen Years = _n-13 if _n<=25
gen Zero = 0 if _n <=25
gen b = 0
gen se = 0
gen u = 0 
gen d = 0

* Run regressions for the likelihood of a complaint investigation before and after citizen complaint
foreach h in 0 1 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_complaint_inv_`h'  p_air_incident, absorb(RN_year t) cluster(RN_id)
replace b = _b[p_air_incident] if Years == `h'
replace se = _se[p_air_incident] if Years == `h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == `h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == `h'
}
foreach h in 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_complaint_inv_neg`h'  p_air_incident, absorb(RN_year t) cluster(RN_id)
replace b = _b[p_air_incident] if Years == -`h'
replace se = _se[p_air_incident] if Years == -`h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == -`h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == -`h'
}

* Save point estimates output
keep if Years != .
keep b u d se Years Zero
export delimited "$point_estimates/Point_Estimates_Figure_3_Panel_C.csv", replace

* Create Figure
graph set window fontface "Times New Roman"
twoway(rarea u d Years, col(gs10) fint(inten30) lwidth(0) lpattern(solid)) ///
    (line b Years, lcolor(black) lpattern(solid) lwidth(medium)) ///
    (line Zero Years, lcolor(gs8)), ///
    xlabel(-12(1)12, nogrid labsize(large)) ///
    ylabel(0(0.1)0.6, labsize(large)) ///
    legend(off) ///
    ytitle("{&Delta} P(Complaint Investigation)", size(large)) ///
    xtitle("Month", size(large)) ///
    graphregion(color(white)) ///
    plotregion(color(white)) ///
    xsize(8.6)
graph export "$figures/Figure_3_Panel_C.pdf", replace


* Panel D: Air Non-Complaint Investigations
use "$processed_data/Air_Panel.dta", clear
drop if never_air_inv == 1

egen t = group(year month)
xtset RN_id t

* Generate difference variables
forv h = 0/12 {
	gen p_air_nocomplaint_inv_`h' = f`h'.p_air_nocomplaint_inv - l1.p_air_nocomplaint_inv
}

forv h = 2/12 {
	gen p_air_nocomplaint_inv_neg`h' = l`h'.p_air_nocomplaint_inv - l1.p_air_nocomplaint_inv
}

egen RN_year = group(RN year)

keep RN_id t RN_year p_air_nocomplaint_inv* p_air_incident

* Initialize storage
cap drop b u d se Years Zero
gen Years = _n-13 if _n<=25
gen Zero = 0 if _n <=25
gen b = 0
gen se = 0
gen u = 0 
gen d = 0

* Run regressions for likelihood of noncomplaint investigation before and after a citizen complaint, given the complaint is from outside TCEQ's jurisdiction.
foreach h in 0 1 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_nocomplaint_inv_`h'  p_air_incident, absorb(RN_year t) cluster(RN_id)
replace b = _b[p_air_incident] if Years == `h'
replace se = _se[p_air_incident] if Years == `h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == `h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == `h'
}
foreach h in 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_nocomplaint_inv_neg`h'  p_air_incident, absorb(RN_year t) cluster(RN_id)
replace b = _b[p_air_incident] if Years == -`h'
replace se = _se[p_air_incident] if Years == -`h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == -`h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == -`h'
}

* Save point estimates output
keep if Years != .
keep b u d se Years Zero
export delimited "$point_estimates/Point_Estimates_Figure_3_Panel_D.csv", replace
graph set window fontface "Times New Roman"
twoway(rarea u d Years, col(gs10) fint(inten30) lwidth(0) lpattern(solid)) ///
    (line b Years, lcolor(gs3) lpattern(solid) lwidth(medium)) ///
    (line Zero Years, lcolor(gs8)), ///
    xlabel(-12(1)12, nogrid labsize(vlarge)) ///
    ylabel(0(0.01)0.06, labsize(vlarge)) ///
    legend(off) ///
    ytitle("{&Delta} P(Non-Complaint Investigation)", size(vlarge)) ///
    xtitle("Month", size(vlarge)) ///
    graphregion(color(white)) ///
    plotregion(color(white)) ///
    xsize(8.6)

* Create Figure
graph export "$figures/Figure_3_Panel_D.pdf", replace